/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package genomemap.worker.helper;

import javautil.collections.ArrayUtil;
import javautil.lang.MathUtil;

/**
 * Class description
 *
 * @version        Enter version here..., 12/11/23
 * @author         Susanta Tewari
 */
public class GHoodHelper {

    /** Field description */
    private static int count = 0;

    /** Field description */
    private static int[] k = null;

    /** Field description */
    private static int[] counter = null;

    /** Field description */
    private static final int[] firstCol = { 1, 1, 0, 0 };;

    /** Field description */
    private static int kCount1 = 0;

    /** Field description */
    private static final int[][] tempVec = {

        { 1, 1, 0, 0 }, { 0, 1, 1, 0 }, { 1, 0, 1, 0 }, { 1, 0, 0, 1 }, { 0, 1, 0, 1 },
        { 0, 0, 1, 1 }
    };;

    /** Field description */
    private static int[] kCount = null;

    public static int getSegments(int loci, int height) {

        if ((loci - 1) % (height - 1) != 0) {

            System.out.println("The height should be adjusted to the loci length");
            System.exit(0);
        }

        return (loci - 1) / (height - 1);
    }

    public static void getActiveXoverFirst(int height, int[][] genotype, MyInteger intObj) {
        ActiveKLoopFirst(0, height, getGenotypeSegment(0, height - 1, genotype), intObj);
    }

    /**
     * Method description
     *
     * @param index description
     * @param height description
     * @param data description
     * @param intObj description
     */
    private static void ActiveKLoopFirst(int index, int height, int[][] data, MyInteger intObj) {

        if (index == 0) {

            count   = 0;
            k       = new int[height - 1];
            counter = new int[height - 2];
        }

        if (index < height - 2) {

            for (counter[index] = 0; counter[index] < 24; counter[index]++) {

                ActiveKLoopFirst(index + 1, height, data, intObj);
            }

            intObj.value = count;
        }

        if (index == height - 2) {

            for (int check1 = 0; check1 < 24; check1++) {

                for (int check2 = 0; check2 < height - 2; check2++) {

                    k[check2] = counter[check2];
                }

                k[height - 2] = check1;

                int[][] r;

                r = getRMatrix(k, firstCol);

                if (kIsWorthy(r, data) == 1) {
                    count++;
                }
            }
        }
    }

    public static void getInfo4First(int height, int[][] genotype, int[][] kArrayFirst,
                                     int[] linkInfoFirst, int[][][] cArrayFirst) {

        int[][] data = getGenotypeSegment(0, height - 1, genotype);

        LoopFirst(0, height, data, kArrayFirst, linkInfoFirst, cArrayFirst);
    }

    /**
     * Method description
     *
     * @param index description
     * @param height description
     * @param data description
     * @param kArrayFirst description
     * @param linkInfoFirst description
     * @param cArrayFirst description
     */
    private static void LoopFirst(int index, int height, int[][] data, int[][] kArrayFirst,
                                  int[] linkInfoFirst, int[][][] cArrayFirst) {

        if (index == 0) {

            kCount1 = 0;
            k       = new int[height - 1];
            counter = new int[height - 2];
        }

        if (index < height - 2) {

            for (counter[index] = 0; counter[index] < 24; counter[index]++) {

                LoopFirst(index + 1, height, data, kArrayFirst, linkInfoFirst, cArrayFirst);
            }
        }

        if (index == height - 2) {

            for (int check1 = 0; check1 < 24; check1++) {

                for (int check2 = 0; check2 < height - 2; check2++) {

                    k[check2] = counter[check2];
                }

                k[height - 2] = check1;

                int[][] r;

                r = getRMatrix(k, firstCol);

                if (kIsWorthy(r, data) == 1) {

                    kArrayFirst[kCount1]   = k.clone();
                    linkInfoFirst[kCount1] = Map6ToInt(r[r.length - 1]);

                    for (int index1 = 0; index1 < data[0].length; index1++) {

                        cArrayFirst[kCount1][index1] = MatchArray(r,
                                ArrayUtil.getCol(data, index1 + 1));
                    }

                    kCount1++;
                }
            }
        }
    }

    public static void getActiveXoverLastArray(int height, int[][] genotype, int[] activeKLast) {

        ActiveKLoopLastArray(0, height,
                             getGenotypeSegment(genotype.length - height, genotype.length - 1,
                                 genotype), activeKLast);
    }

    public static void getInfo4Last(int height, int[][] genotype, int[][][] kArrayLast,
                                    int[][][][] cArrayLast) {

        int[][] data = getGenotypeSegment(genotype.length - height, genotype.length - 1, genotype);

        LoopLast(0, height, data, kArrayLast, cArrayLast);
    }

    /**
     * Method description
     *
     * @param index description
     * @param height description
     * @param data description
     * @param activeKLast description
     */
    private static void ActiveKLoopLastArray(int index, int height, int[][] data,
            int[] activeKLast) {

        if (index == 0) {

            k       = new int[height - 1];
            counter = new int[height - 2];
        }

        if (index < height - 2) {

            for (counter[index] = 0; counter[index] < 24; counter[index]++) {

                ActiveKLoopLastArray(index + 1, height, data, activeKLast);
            }
        }

        if (index == height - 2) {

            for (int check1 = 0; check1 < 24; check1++) {

                for (int check2 = 0; check2 < height - 2; check2++) {

                    k[check2] = counter[check2];
                }

                k[height - 2] = check1;

                int[][] r;

                for (int i = 0; i < 6; i++) {

                    r = getRMatrix(k, tempVec[i]);    // r is  segment-1 by 4

                    if (kIsWorthy(r, data) == 1) {
                        activeKLast[i]++;
                    }
                }
            }
        }
    }

    /**
     * Method description
     *
     * @param index description
     * @param height description
     * @param data description
     * @param kArrayLast description
     * @param cArrayLast description
     */
    private static void LoopLast(int index, int height, int[][] data, int[][][] kArrayLast,
                                 int[][][][] cArrayLast) {

        if (index == 0) {

            kCount  = new int[6];
            k       = new int[height - 1];
            counter = new int[height - 2];
        }

        if (index < height - 2) {

            for (counter[index] = 0; counter[index] < 24; counter[index]++) {

                LoopLast(index + 1, height, data, kArrayLast, cArrayLast);
            }
        }

        if (index == height - 2) {

            for (int check1 = 0; check1 < 24; check1++) {

                for (int check2 = 0; check2 < height - 2; check2++) {

                    k[check2] = counter[check2];
                }

                k[height - 2] = check1;

                int[][] r;

                for (int i = 0; i < 6; i++) {

                    r = getRMatrix(k, tempVec[i]);    // r is  segment-1 by 4

                    if (kIsWorthy(r, data) == 1) {

                        kArrayLast[Map6ToInt(tempVec[i])][kCount[i]] = k.clone();

                        for (int index1 = 0; index1 < data[0].length; index1++) {

                            cArrayLast[Map6ToInt(tempVec[i])][kCount[i]][index1] = MatchArray(r,
                                    ArrayUtil.getCol(data, index1 + 1));
                        }

                        kCount[i]++;
                    }
                }
            }
        }
    }

    public static void ActiveKLoopTemp(int index, int height, int indexS, int[] activeKTemp,
                                       int[][] data) {

        if (index == 0) {

            k       = new int[height - 1];
            counter = new int[height - 2];
        }

        if (index < height - 2) {

            for (counter[index] = 0; counter[index] < 24; counter[index]++) {

                ActiveKLoopTemp(index + 1, height, indexS, activeKTemp, data);
            }
        }

        if (index == height - 2) {

            for (int check1 = 0; check1 < 24; check1++) {

                for (int check2 = 0; check2 < height - 2; check2++) {

                    k[check2] = counter[check2];
                }

                k[height - 2] = check1;

                int[][] r;

                for (int i = 0; i < 6; i++) {

                    r = getRMatrix(k, tempVec[i]);    // r is  segment-1 by 4

                    if (kIsWorthy(r, data) == 1) {
                        activeKTemp[i]++;
                    }
                }
            }
        }
    }

    public static void LoopTemp(int index, int height, int indexS, int[][] data,
                                int[][][] kArrayTemp, int[][] linkInfoTemp,
                                int[][][][] cArrayTemp) {

        if (index == 0) {

            kCount  = new int[6];
            k       = new int[height - 1];
            counter = new int[height - 2];
        }

        if (index < height - 2) {

            for (counter[index] = 0; counter[index] < 24; counter[index]++) {

                LoopTemp(index + 1, height, indexS, data, kArrayTemp, linkInfoTemp, cArrayTemp);
            }
        }

        if (index == height - 2) {

            for (int check1 = 0; check1 < 24; check1++) {

                for (int check2 = 0; check2 < height - 2; check2++) {

                    k[check2] = counter[check2];
                }

                k[height - 2] = check1;

                int[][] r;

                /*
                 * -----------------------------------------------------
                 *      The following code generates cArray for all possible
                 *      16 columns in the end.
                 * ------------------------------------------------
                 */
                for (int i = 0; i < 6; i++) {

                    r = getRMatrix(k, tempVec[i]);    // r is  segment-1 by 4

                    if (kIsWorthy(r, data) == 1) {

                        kArrayTemp[Map6ToInt(tempVec[i])][kCount[i]]   = k.clone();
                        linkInfoTemp[Map6ToInt(tempVec[i])][kCount[i]] = Map6ToInt(r[r.length - 1]);

                        for (int index1 = 0; index1 < data[0].length; index1++) {

                            cArrayTemp[Map6ToInt(tempVec[i])][kCount[i]][index1] = MatchArray(r,
                                    ArrayUtil.getCol(data, index1 + 1));
                        }

                        kCount[i]++;
                    }
                }
            }
        }
    }

    public static int[][] getGenotypeSegment(int startRowIndex, int endRowIndex, int[][] genotype) {

        int[][] result = new int[endRowIndex - startRowIndex + 1][genotype[0].length];
        int count1     = 0;

        for (int i = startRowIndex; i <= endRowIndex; i++) {

            for (int j = 0; j < genotype[0].length; j++) {

                result[count1][j] = genotype[i][j];
            }

            count1++;
        }

        return result;
    }

    /**
     * Method description
     *
     * @param i description
     *
     * @return description
     */
    private static int CellSpecial(int i) {

        int result = 0;

        if (i == 0) {
            result = 0;
        } else if (((i >= 1) && (i <= 5)) || (i == 10) || (i == 15) || (i == 20)) {
            result = 1;
        } else {
            result = 2;
        }

        return result;
    }

    /**
     * Method description
     *
     * @param loci description
     * @param height description
     * @param distinctGenotypes description
     * @param genotype description
     * @param genotypeFreq description
     * @param cProbOld description
     * @param activeKFirst description
     * @param kArrayFirst description
     * @param linkInfoFirst description
     * @param cArrayFirst description
     * @param activeKLast description
     * @param kArrayLast description
     * @param cArrayLast description
     * @param spores description
     * @param tempSumIndex description
     * @param probFOld description
     * @param pCount description
     * @param postCount description
     */
    private static void SpiralUp(int loci, int height, int distinctGenotypes, int[][] genotype,
                                 int[] genotypeFreq, double[] cProbOld, int activeKFirst,
                                 int[][] kArrayFirst, int[] linkInfoFirst, int[][][] cArrayFirst,
                                 int[] activeKLast, int[][][] kArrayLast, int[][][][] cArrayLast,
                                 int[][] spores, int[][][][] tempSumIndex, double[] probFOld,
                                 double[][][] pCount, double[][] postCount) {

        double[][][] FirstTempSum        = new double[6][distinctGenotypes][11];
        double[][][][][] FirstTempPCount = new double[6][distinctGenotypes][11][3][loci - 1];

        /*
         * ----------------------------------
         *  Create the first tempSum.
         * ---------------------------------
         */
        double[] lastCProb = getCProb(cProbOld, 0, height - 1);

        for (int index0 = 0; index0 < 6; index0++) {

            for (int index1 = 0; index1 < distinctGenotypes; index1++) {

                for (int index3 = 0; index3 < 11; index3++) {

                    for (int index2 = 0; index2 < activeKLast[index0]; index2++) {

                        double tempDouble = kProb(lastCProb, kArrayLast[index0][index2]);
                        double condProb   = tempSumIndex[index0][index1][index3][index2] / 4.0;

                        FirstTempSum[index0][index1][index3] += condProb * tempDouble;

                        for (int index = loci - height; index < loci - 1; index++) {

                            FirstTempPCount[index0][index1][index3][CellSpecial(kArrayLast[index0][index2][index - loci + height])][index] +=
                                condProb * tempDouble;
                        }
                    }
                }
            }
        }

        double[][][] temp1TempSum        = FirstTempSum;
        double[][][][][] temp1TempPCount = FirstTempPCount;
        double[][][] temp2TempSum        = new double[6][distinctGenotypes][11];
        double[][][][][] temp2TempPCount = new double[6][distinctGenotypes][11][3][loci - 1];

        /*
         * ----------------------------
         * Intermediate steps.
         * ---------------------------
         */
        int startRowIndex = 0;
        int endRowIndex   = loci - height;    // note that there is a common gene in each interval
        int startPos      = 0;
        int endPos        = loci - height;
        int segments      = getSegments(loci, height);

        for (int indexS = 0; indexS < segments - 2; indexS++) {

            startPos = endPos - height + 2;

            double[] firstCProb = ChopAtoB(cProbOld, startPos, endPos);


            // setting up the sth temporary segment
            int[] activeKTemp = new int[6];


            // Get the relevant data
            startRowIndex = endRowIndex - height + 1;    // endRowIndex is changed in the bottom.

            int[][] data = getGenotypeSegment(startRowIndex, endRowIndex, genotype);


            // Calculate active crossovers. Stored in activeKTemp.
            ActiveKLoopTemp(0, height, indexS, activeKTemp, data);    // note this uses data.


            // allocate cells.
            int[][][][] cArrayTemp = new int[6][][][];
            int[][] linkInfoTemp   = new int[6][];
            int[][][] kArrayTemp   = new int[6][][];

            for (int i = 0; i < 6; i++) {

                cArrayTemp[i]   = new int[activeKTemp[i]][distinctGenotypes][4];
                linkInfoTemp[i] = new int[activeKTemp[i]];
                kArrayTemp[i]   = new int[activeKTemp[i]][height - 1];
            }


            // run the loop to set kArrayTemp,cArrayTemp,linkInfoTemp.
            LoopTemp(0, height, indexS, data, kArrayTemp, linkInfoTemp, cArrayTemp);    // note this uses data


            // linking at work
            for (int index0 = 0; index0 < 6; index0++) {

                for (int index1 = 0; index1 < distinctGenotypes; index1++) {

                    for (int index2 = 0; index2 < 11; index2++) {

                        for (int indexK = 0; indexK < activeKTemp[index0]; indexK++) {

                            int[] k     = kArrayTemp[index0][indexK].clone();
                            double prob = kProb(firstCProb, k);
                            int index01 = linkInfoTemp[index0][indexK];
                            int spike =
                                getSporeMatch(UpdateSpore(cArrayTemp[index0][indexK][index1],
                                                          spores[index2]));

                            temp2TempSum[index0][index1][index2] +=
                                temp1TempSum[index01][index1][spike] * prob;

                            for (int index = startPos - 1; index < endPos; index++) {

                                temp2TempPCount[index0][index1][index2][CellSpecial(k[index - startPos + 1])][index] +=
                                    temp1TempSum[index01][index1][spike] * prob;
                            }

                            for (int i = 0; i < 3; i++) {

                                for (int index = endPos; index < loci - 1; index++) {

                                    temp2TempPCount[index0][index1][index2][i][index] +=
                                        temp1TempPCount[index01][index1][spike][i][index] * prob;
                                }
                            }
                        }
                    }
                }
            }


            // change temp1 to temp2
            temp1TempSum    = temp2TempSum.clone();
            temp1TempPCount = temp2TempPCount.clone();


            // set temp2 to zeros for security.
            temp2TempSum    = new double[6][distinctGenotypes][11];
            temp2TempPCount = new double[6][distinctGenotypes][11][3][loci - 1];


            // System.gc(); // call grabage collection
            endPos      = startPos - 1;
            endRowIndex = startRowIndex;
        }

        /*
         * ----------------------------
         * Final Loop
         * ---------------------------
         */
        double[] firstCProb = getCProb(cProbOld, 1, height - 1);
        int total_obs       = MathUtil.Sum(genotypeFreq);

        for (int index1 = 0; index1 < activeKFirst; index1++) {

            int[] k     = kArrayFirst[index1].clone();
            double prob = kProb(firstCProb, k);
            int[] spike = new int[distinctGenotypes];


            // spike is updated spore before the last segment
            for (int i = 0; i < distinctGenotypes; i++) {

                spike[i] = getSporeMatch(cArrayFirst[index1][i]);
            }


            // probFOld and pCount are global variables.
            int index0 = linkInfoFirst[index1];

            for (int index3 = 0; index3 < distinctGenotypes; index3++) {

                probFOld[index3] += temp1TempSum[index0][index3][spike[index3]] * prob;

                for (int index = 0; index < height - 1; index++) {

                    pCount[CellSpecial(k[index])][index][index3] +=
                        temp1TempSum[index0][index3][spike[index3]] * prob * genotypeFreq[index3]
                        / total_obs;
                }

                for (int i = 0; i < 3; i++) {

                    for (int index = height - 1; index < loci - 1; index++) {

                        pCount[i][index][index3] +=
                            temp1TempPCount[index0][index3][spike[index3]][i][index] * prob
                            * genotypeFreq[index3] / total_obs;
                    }
                }
            }
        }
    }

    public static double[] getCProbNew(int loci, int height, int distinctGenotypes,
                                       int[][] genotype, int[] genotypeFreq, double[] cProbOld,
                                       int activeKFirst, int[][] kArrayFirst, int[] linkInfoFirst,
                                       int[][][] cArrayFirst, int[] activeKLast,
                                       int[][][] kArrayLast, int[][][][] cArrayLast,
                                       int[][] spores, int[][][][] tempSumIndex,
                                       double[] probFOld) {

        double[][][] pCount  = new double[3][loci - 1][distinctGenotypes];
        double[][] postCount = new double[3][loci - 1];

        SpiralUp(loci, height, distinctGenotypes, genotype, genotypeFreq, cProbOld, activeKFirst,
                 kArrayFirst, linkInfoFirst, cArrayFirst, activeKLast, kArrayLast, cArrayLast,
                 spores, tempSumIndex, probFOld, pCount, postCount);

        return calculateCProb(loci, distinctGenotypes, postCount, pCount, probFOld);
    }

    /**
     * Method description
     *
     * @param loci description
     * @param distinctGenotypes description
     * @param postCount description
     * @param pCount description
     * @param probFOld description
     *
     * @return description
     */
    private static double[] calculateCProb(int loci, int distinctGenotypes, double[][] postCount,
            double[][][] pCount, double[] probFOld) {

        double[] cProbNew = new double[loci - 1];

        for (int index1 = 0; index1 < 3; index1++) {

            for (int index2 = 0; index2 < loci - 1; index2++) {

                for (int index3 = 0; index3 < distinctGenotypes; index3++) {

                    pCount[index1][index2][index3] = pCount[index1][index2][index3]
                                                     / probFOld[index3];
                }

                postCount[index1][index2] = MathUtil.Sum(pCount[index1][index2]);
            }
        }

        for (int index1 = 0; index1 < loci - 1; index1++) {

            cProbNew[index1] = (postCount[1][index1] + 2 * postCount[2][index1])
                               / (2.0 * MathUtil.Sum(ArrayUtil.getCol(postCount, index1 + 1)));
        }

        return cProbNew;
    }

    public static int[][][][] getTempSumIndex(int distinctGenotypes, int[][] spores,
            int[] activeKLast, int[][][][] cArrayLast) {

        int[][][][] tempSumIndex = new int[6][][][];

        for (int i = 0; i < 6; i++) {

            tempSumIndex[i] = new int[distinctGenotypes][11][activeKLast[i]];
        }

        for (int index0 = 0; index0 < 6; index0++) {

            for (int index1 = 0; index1 < distinctGenotypes; index1++) {

                for (int index2 = 0; index2 < 11; index2++) {

                    for (int index3 = 0; index3 < activeKLast[index0]; index3++) {

                        tempSumIndex[index0][index1][index2][index3] =
                            MathUtil.Sum(UpdateSpore(cArrayLast[index0][index3][index1],
                                                     spores[index2]));
                    }
                }
            }
        }

        return tempSumIndex;
    }

    public static double kProbDeriOff(double[] cProbOld, int[] k, int i, int j) {

        double product = 1;

        for (int index = 0; index < k.length; index++) {

            if ((index != (i - 1)) && (index != (j - 1))) {
                product *= lociProb(k[index], cProbOld[index]);
            } else {

                if (index == i - 1) {
                    product = product * probD(cProbOld, k, i);
                } else {
                    product = product * probD(cProbOld, k, j);
                }
            }
        }

        return product;
    }

    /**
     * Method description
     *
     * @param cProbOld description
     * @param k description
     * @param i description
     *
     * @return description
     */
    private static double probD(double[] cProbOld, int[] k, int i) {

        double result = 0.0;

        if (k[i - 1] == 0) {
            result = 2 * (1 - cProbOld[i - 1]) * (-1);
        } else if (((k[i - 1] >= 1) && (k[i - 1] <= 5)) || (k[i - 1] == 10) || (k[i - 1] == 15)
                   || (k[i - 1] == 20)) {
            result = 0.25 - cProbOld[i - 1] / 2.0;
        } else {
            result = cProbOld[i - 1] / 8.0;
        }

        return result;
    }

    public static double kProbDeriDiag(double[] cProbOld, int[] k, int i) {

        double product = 1;

        for (int index = 0; index < k.length; index++) {

            if (index != (i - 1)) {
                product *= lociProb(k[index], cProbOld[index]);
            } else {
                product = product * probD2(cProbOld, k, i);
            }
        }

        return product;
    }

    /**
     * Method description
     *
     * @param cProbOld description
     * @param k description
     * @param i description
     *
     * @return description
     */
    private static double probD2(double[] cProbOld, int[] k, int i) {

        double result = 0.0;

        if (k[i - 1] == 0) {
            result = 2;
        } else if (((k[i - 1] >= 1) && (k[i - 1] <= 5)) || (k[i - 1] == 10) || (k[i - 1] == 15)
                   || (k[i - 1] == 20)) {
            result = -0.5;
        } else {
            result = 1.0 / 8.0;
        }

        return result;
    }

    public static double kProbSingle(double[] cProbOld, int[] k, int i) {

        double product = 1;

        for (int index = 0; index < k.length; index++) {

            if (index != (i - 1)) {
                product *= lociProb(k[index], cProbOld[index]);
            } else {
                product = product * probD(cProbOld, k, i);
            }
        }

        return product;
    }

    public static int[] ChopAtoB(int[] array, int pos1, int pos2) {

        /*
         * -------------------------------------------------
         * chops A to B including themselves.
         * Note that, pos1 and pos2 are positions not indices
         * and output is from the small pos to the bigger one.
         * --------------------------------------------------
         */
        int[] result = { 1 };

        if (pos1 > pos2) {

            result = new int[pos1 - pos2 + 1];

            for (int index = 0; index < result.length; index++) {

                result[index] = array[pos2 - 1 + index];
            }
        }

        if (pos1 < pos2) {

            result = new int[pos2 - pos1 + 1];

            for (int index = 0; index < result.length; index++) {

                result[index] = array[pos1 - 1 + index];
            }
        }

        if (pos1 == pos2) {
            result[0] = array[pos1 - 1];
        }

        return result;
    }

    public static double[] ChopAtoB(double[] array, int pos1, int pos2) {

        /*
         * -------------------------------------------------
         * chops A to B including themselves.
         * Note that, pos1 and pos2 are positions not indices
         * and output is from the small pos to the bigger one.
         * --------------------------------------------------
         */
        double[] result = { 1 };

        if (pos1 > pos2) {

            result = new double[pos1 - pos2 + 1];

            for (int index = 0; index < result.length; index++) {

                result[index] = array[pos2 - 1 + index];
            }
        }

        if (pos1 < pos2) {

            result = new double[pos2 - pos1 + 1];

            for (int index = 0; index < result.length; index++) {

                result[index] = array[pos1 - 1 + index];
            }
        }

        if (pos1 == pos2) {
            result[0] = array[pos1 - 1];
        }

        return result;
    }

    /**
     * Method description
     *
     * @param cProb description
     * @param index description
     * @param length description
     *
     * @return description
     */
    private static double[] getCProb(double[] cProb, int index, int length) {

        double[] result = new double[length];

        if (index == 1) {    // front slice

            for (int index1 = 0; index1 < length; index1++) {

                result[index1] = cProb[index1];
            }
        }

        if (index == 0) {    // back slice

            int total = cProb.length;

            for (int index2 = 0; index2 < length; index2++) {

                result[index2] = cProb[total - length + index2];
            }
        }

        return result;
    }

    public static int[][] getSpores() {

        int[][] result = {

            { 0, 0, 0, 0 }, { 1, 0, 0, 0 }, { 0, 1, 0, 0 }, { 0, 0, 1, 0 }, { 0, 0, 0, 1 },
            { 1, 1, 0, 0 }, { 1, 0, 1, 0 }, { 1, 0, 0, 1 }, { 0, 1, 1, 0 }, { 0, 1, 0, 1 },
            { 0, 0, 1, 1 }
        };

        return result;
    }

    public static int getSporeMatch(int[] spore) {

        int result   = 13;    // for guarding against error.
        int[][] temp = getSpores();

        for (int index = 0; index < temp.length; index++) {

            if (isMatch(temp[index], spore) == 1) {

                result = index;

                break;
            }
        }

        return result;
    }

    public static int[] UpdateSpore(int[] array1, int[] array2) {


        // in the application array1 and array2 both have 4 columns.
        // rows denote all possible genotypes.
        if (!((array1.length == 4) && (array2.length == 4))) {

            System.out.println("in UpdateSpore colns don`t have length 4");
            System.exit(0);
        }

        int[] result = new int[4];

        for (int index1 = 0; index1 < array1.length; index1++) {

            if ((array1[index1] == 1) && (array2[index1] == 1)) {
                result[index1] = 1;
            } else {
                result[index1] = 0;
            }
        }

        return result;
    }

    public static double kProb(double[] c, int[] order) {

        double product = 1;

        for (int index = 0; index < order.length; index++) {

            product *= lociProb(order[index], c[index]);
        }

        return product;
    }

    /**
     * Method description
     *
     * @param array description
     *
     * @return description
     */
    private static int Map6ToInt(int[] array) {

        int result = 10;    // this is set to check error.As this would give index out


        // of bound error in the calling place.
        int[] vec0 = { 1, 1, 0, 0 };
        int[] vec1 = { 0, 1, 1, 0 };
        int[] vec2 = { 1, 0, 1, 0 };
        int[] vec3 = { 1, 0, 0, 1 };
        int[] vec4 = { 0, 1, 0, 1 };
        int[] vec5 = { 0, 0, 1, 1 };

        if (isMatch(array, vec0) == 1) {
            result = 0;
        } else if (isMatch(array, vec1) == 1) {
            result = 1;
        } else if (isMatch(array, vec2) == 1) {
            result = 2;
        } else if (isMatch(array, vec3) == 1) {
            result = 3;
        } else if (isMatch(array, vec4) == 1) {
            result = 4;
        } else if (isMatch(array, vec5) == 1) {
            result = 5;
        }

        return result;
    }

    public static double Convergence(double[] array1, double[] array2, double limit) {

        /*
         *  double bigone,smallone;
         * bigone=Math.abs(array1[0]-array2[0]);
         * for(int index=1;index<array1.length;index++){
         *    smallone=Math.abs(array1[index]-array2[index]);
         *    if(smallone > bigone){
         *        bigone=smallone;
         *    }
         * }
         */

        /* flag=1 means convergence achieved and not otherwise. */
        int flag = 1;

        for (int index = 0; index < array1.length; index++) {

            if (Math.abs(array1[index] - array2[index]) > limit) {
                flag = 0;
            }
        }

        return flag;
    }

    /**
     * Method description
     *
     * @param array1 description
     * @param array2 description
     *
     * @return description
     */
    private static int kIsWorthy(int[][] array1, int[][] array2) {

        /* result =1 means atleast one match is there and 0 means none */
        int result = 0;

        for (int index = 0; index < array2[0].length; index++) {

            if (countMatch(array1, ArrayUtil.getCol(array2, index + 1)) > 0) {

                result = 1;

                break;
            }
        }

        return result;
    }

    /**
     * Method description
     *
     * @param array description
     * @param vector description
     *
     * @return description
     */
    private static int countMatch(int[][] array, int[] vector) {

        int count = 0;

        for (int index = 0; index < array[0].length; index++) {

            count += isMatch(ArrayUtil.getCol(array, index + 1), vector);
        }

        return count;
    }

    /**
     * Method description
     *
     * @param array description
     * @param vector description
     *
     * @return description
     */
    private static int[] MatchArray(int[][] array, int[] vector) {


        // Note that array is somesize by 4, so there`re only 4 colns.
        int[] result = new int[4];

        for (int index = 0; index < array[0].length; index++) {

            result[index] = isMatch(ArrayUtil.getCol(array, index + 1), vector);
        }

        return result;
    }

    /**
     * Method description
     *
     * @param array1 description
     * @param array2 description
     *
     * @return description
     */
    private static int isMatch(int[] array1, int[] array2) {

        /* returning result=1 means a match and 0 a mismatch */
        int result = 1;

        for (int index = 0; index < array1.length; index++) {

            if (array1[index] != array2[index]) {

                result = 0;

                break;
            }
        }

        return result;
    }

    /**
     * Method description
     *
     * @param array description
     * @param col description
     *
     * @return description
     */
    private static int[][] getRMatrix(int[] array, int[] col) {


        // array is k: and col is the first column
        int[][] temp = new int[array.length + 1][4];

        temp[0][0] = col[0];
        temp[0][1] = col[1];
        temp[0][2] = col[2];
        temp[0][3] = col[3];


        // getTPattern is creating the next tetrad pattern
        for (int index = 1; index <= array.length; index++) {

            switch (array[index - 1]) {

                case 0 :
                    temp[index] = getTPattern(temp[index - 1], 0, 0);

                    break;

                case 1 :
                    temp[index] = getTPattern(temp[index - 1], 0, 1);

                    break;

                case 2 :
                    temp[index] = getTPattern(temp[index - 1], 0, 2);

                    break;

                case 3 :
                    temp[index] = getTPattern(temp[index - 1], 0, 3);

                    break;

                case 4 :
                    temp[index] = getTPattern(temp[index - 1], 0, 4);

                    break;

                case 5 :
                    temp[index] = getTPattern(temp[index - 1], 1, 0);

                    break;

                case 6 :
                    temp[index] = getTPattern(temp[index - 1], 1, 1);

                    break;

                case 7 :
                    temp[index] = getTPattern(temp[index - 1], 1, 2);

                    break;

                case 8 :
                    temp[index] = getTPattern(temp[index - 1], 1, 3);

                    break;

                case 9 :
                    temp[index] = getTPattern(temp[index - 1], 1, 4);

                    break;

                case 10 :
                    temp[index] = getTPattern(temp[index - 1], 2, 0);

                    break;

                case 11 :
                    temp[index] = getTPattern(temp[index - 1], 2, 1);

                    break;

                case 12 :
                    temp[index] = getTPattern(temp[index - 1], 2, 2);

                    break;

                case 13 :
                    temp[index] = getTPattern(temp[index - 1], 2, 3);

                    break;

                case 14 :
                    temp[index] = getTPattern(temp[index - 1], 2, 4);

                    break;

                case 15 :
                    temp[index] = getTPattern(temp[index - 1], 3, 0);

                    break;

                case 16 :
                    temp[index] = getTPattern(temp[index - 1], 3, 1);

                    break;

                case 17 :
                    temp[index] = getTPattern(temp[index - 1], 3, 2);

                    break;

                case 18 :
                    temp[index] = getTPattern(temp[index - 1], 3, 3);

                    break;

                case 19 :
                    temp[index] = getTPattern(temp[index - 1], 3, 4);

                    break;

                case 20 :
                    temp[index] = getTPattern(temp[index - 1], 4, 0);

                    break;

                case 21 :
                    temp[index] = getTPattern(temp[index - 1], 4, 1);

                    break;

                case 22 :
                    temp[index] = getTPattern(temp[index - 1], 4, 2);

                    break;

                case 23 :
                    temp[index] = getTPattern(temp[index - 1], 4, 3);

                    break;

                case 24 :
                    temp[index] = getTPattern(temp[index - 1], 4, 4);

                    break;
            }
        }

        return temp;
    }

    public static double lociProb(int k, double c) {

        double result = 1;

        switch (k) {

            case 0 :
                result = SetProb(c, 0, 0);

                break;

            case 1 :
                result = SetProb(c, 0, 1);

                break;

            case 2 :
                result = SetProb(c, 0, 2);

                break;

            case 3 :
                result = SetProb(c, 0, 3);

                break;

            case 4 :
                result = SetProb(c, 0, 4);

                break;

            case 5 :
                result = SetProb(c, 1, 0);

                break;

            case 6 :
                result = SetProb(c, 1, 1);

                break;

            case 7 :
                result = SetProb(c, 1, 2);

                break;

            case 8 :
                result = SetProb(c, 1, 3);

                break;

            case 9 :
                result = SetProb(c, 1, 4);

                break;

            case 10 :
                result = SetProb(c, 2, 0);

                break;

            case 11 :
                result = SetProb(c, 2, 1);

                break;

            case 12 :
                result = SetProb(c, 2, 2);

                break;

            case 13 :
                result = SetProb(c, 2, 3);

                break;

            case 14 :
                result = SetProb(c, 2, 4);

                break;

            case 15 :
                result = SetProb(c, 3, 0);

                break;

            case 16 :
                result = SetProb(c, 3, 1);

                break;

            case 17 :
                result = SetProb(c, 3, 2);

                break;

            case 18 :
                result = SetProb(c, 3, 3);

                break;

            case 19 :
                result = SetProb(c, 3, 4);

                break;

            case 20 :
                result = SetProb(c, 4, 0);

                break;

            case 21 :
                result = SetProb(c, 4, 1);

                break;

            case 22 :
                result = SetProb(c, 4, 2);

                break;

            case 23 :
                result = SetProb(c, 4, 3);

                break;

            case 24 :
                result = SetProb(c, 4, 4);

                break;
        }

        return result;
    }

    /**
     * Method description
     *
     * @param c description
     * @param i description
     * @param j description
     *
     * @return description
     */
    private static double SetProb(double c, int i, int j) {

        double prod = 1.0;


        // for ith part
        if (i == 0) {
            prod *= (1 - c);
        } else {
            prod *= c / 4.0;
        }


        // for jth part
        if (j == 0) {
            prod *= (1 - c);
        } else {
            prod *= c / 4.0;
        }

        return prod;
    }

    /**
     * Method description
     *
     * @param array description
     * @param i description
     * @param j description
     *
     * @return description
     */
    private static int[] getTPattern(int[] array, int i, int j) {

        int[] step1 = new int[4];    // hard wired to detct problems ..this will always be size 4
        int[] step2 = new int[4];


        // this switch creates step1
        switch (i) {


            // no change
            case 0 :
                step1 = array.clone();

                break;

            case 1 :


                // change 1st and 3rd strand
                step1[0] = array[2];
                step1[2] = array[0];


                // rest are unaffected
                step1[1] = array[1];
                step1[3] = array[3];

                break;

            case 2 :


                // change 2nd and 3rd strand
                step1[1] = array[2];
                step1[2] = array[1];


                // rest are unaffected
                step1[0] = array[0];
                step1[3] = array[3];

                break;

            case 3 :


                // change 2nd and 4th strand
                step1[1] = array[3];
                step1[3] = array[1];


                // rest are unaffected
                step1[0] = array[0];
                step1[2] = array[2];

                break;

            case 4 :


                // change 1st and 4th strand
                step1[0] = array[3];
                step1[3] = array[0];


                // rest are unaffected
                step1[1] = array[1];
                step1[2] = array[2];

                break;
        }


        // this switch creates step2
        switch (j) {


            // no change
            case 0 :
                step2 = step1.clone();

                break;

            case 1 :


                // change 1st and 3rd strand
                step2[0] = step1[2];
                step2[2] = step1[0];


                // rest are unaffected
                step2[1] = step1[1];
                step2[3] = step1[3];

                break;

            case 2 :


                // change 2nd and 3rd strand
                step2[1] = step1[2];
                step2[2] = step1[1];


                // rest are unaffected
                step2[0] = step1[0];
                step2[3] = step1[3];

                break;

            case 3 :


                // change 2nd and 4th strand
                step2[1] = step1[3];
                step2[3] = step1[1];


                // rest are unaffected
                step2[0] = step1[0];
                step2[2] = step1[2];

                break;

            case 4 :


                // change 1st and 4th strand
                step2[0] = step1[3];
                step2[3] = step1[0];


                // rest are unaffected
                step2[1] = step1[1];
                step2[2] = step1[2];

                break;
        }

        return step2;
    }

    public static int[] Concat1D(int[] array1, int[] array2) {

        int[] result = new int[array1.length + array2.length];

        for (int index = 0; index < array1.length; index++) {

            result[index] = array1[index];
        }

        for (int index = 0; index < array2.length; index++) {

            result[array1.length + index] = array2[index];
        }

        return result;
    }
}
